1 Introduction

In this script, we go through all the pre-registered proposed analyses. As a reminder, the main research questions where as follows:

  1. To what extent does infants#’ preference for IDS as measured in a laboratory setting predict their vocabulary at 18 and 24 months?
  2. Does the relation between IDS preference and vocabulary size change over development?
  3. Are there systematic differences in the strength of this relationship across the language communities in our sample?

We first present the main “sample theory based” analyses (also known as frequentist), separately on the North American and UK samples in parallel to answer our first two research questions, then on the total dataset to answer our third research question. We then provide additional Bayesian statistics where a null effect was found, as specified in the pre-registration.

# Library imports, general settings ==============
library(tidyverse); library(egg)
library(knitr)
library(lme4); library(lmerTest); library(simr)
# As in our discussion with Mike, we will use lmerTest for calculating p values
# in the mixed-effects models (noted as a deviation) 
library(brms); library(rstan)
library(future)
plan(multicore, workers = 8) # Swap to multiprocess if running from Rstudio
library(lattice)
library(effects)
library(sjPlot)
library(robustlmm)
library(car)

# Load model comparison functions
source("helper/lrtests.R")

# Deal with package priority issues
select <- dplyr::select

theme_set(theme_bw(base_size = 10))
options("future" = T)
#knitr::opts_chunk$set(cache = TRUE)

print(sessionInfo()) #listing all info about R and packages info
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] furrr_0.2.1          future.apply_1.7.0   bridgesampling_1.0-0
##  [4] car_3.0-10           robustlmm_2.4-4      sjPlot_2.8.7        
##  [7] effects_4.2-0        carData_3.0-4        lattice_0.20-41     
## [10] future_1.21.0        rstan_2.21.2         StanHeaders_2.21.0-7
## [13] brms_2.14.4          Rcpp_1.0.5           simr_1.0.5          
## [16] lmerTest_3.1-3       lme4_1.1-26          Matrix_1.3-2        
## [19] knitr_1.30           egg_0.4.5            gridExtra_2.3       
## [22] forcats_0.5.0        stringr_1.4.0        dplyr_1.0.2         
## [25] purrr_0.3.4          readr_1.4.0          tidyr_1.1.2         
## [28] tibble_3.0.4         ggplot2_3.3.3        tidyverse_1.3.0     
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.1.0    htmlwidgets_1.5.3   grid_4.0.3         
##   [4] munsell_0.5.0       codetools_0.2-18    effectsize_0.4.1   
##   [7] statmod_1.4.35      DT_0.17             miniUI_0.1.1.1     
##  [10] withr_2.3.0         Brobdingnag_1.2-6   colorspace_2.0-0   
##  [13] rstudioapi_0.13     stats4_4.0.3        robustbase_0.93-7  
##  [16] bayesplot_1.8.0     listenv_0.8.0       emmeans_1.5.3      
##  [19] RLRsim_3.1-6        coda_0.19-4         parallelly_1.23.0  
##  [22] vctrs_0.3.6         generics_0.1.0      TH.data_1.0-10     
##  [25] xfun_0.20           R6_2.5.0            markdown_1.1       
##  [28] gamm4_0.2-6         projpred_2.0.2      assertthat_0.2.1   
##  [31] promises_1.1.1      scales_1.1.1        multcomp_1.4-15    
##  [34] nnet_7.3-14         gtable_0.3.0        globals_0.14.0     
##  [37] processx_3.4.5      sandwich_3.0-0      rlang_0.4.10       
##  [40] splines_4.0.3       broom_0.7.3         inline_0.3.17      
##  [43] yaml_2.2.1          reshape2_1.4.4      abind_1.4-5        
##  [46] modelr_0.1.8        threejs_0.3.3       crosstalk_1.1.1    
##  [49] backports_1.2.1     httpuv_1.5.5        rsconnect_0.8.16   
##  [52] tools_4.0.3         ellipsis_0.3.1      ggridges_0.5.3     
##  [55] plyr_1.8.6          base64enc_0.1-3     ps_1.5.0           
##  [58] prettyunits_1.1.1   zoo_1.8-8           haven_2.3.1        
##  [61] fs_1.5.0            survey_4.0          magrittr_2.0.1     
##  [64] data.table_1.13.6   openxlsx_4.2.3      colourpicker_1.1.0 
##  [67] reprex_0.3.0        mvtnorm_1.1-1       sjmisc_2.8.6       
##  [70] matrixStats_0.57.0  hms_1.0.0           shinyjs_2.0.0      
##  [73] mime_0.9            evaluate_0.14       xtable_1.8-4       
##  [76] shinystan_2.5.0     pbkrtest_0.5-0.1    rio_0.5.16         
##  [79] sjstats_0.18.1      readxl_1.3.1        ggeffects_1.0.1    
##  [82] rstantools_2.1.1    compiler_4.0.3      V8_3.4.0           
##  [85] crayon_1.3.4        minqa_1.2.4         htmltools_0.5.1    
##  [88] mgcv_1.8-33         later_1.1.0.1       RcppParallel_5.0.2 
##  [91] lubridate_1.7.9.2   DBI_1.1.0           sjlabelled_1.1.7   
##  [94] dbplyr_1.4.2        MASS_7.3-53         boot_1.3-25        
##  [97] cli_2.2.0           mitools_2.4         parallel_4.0.3     
## [100] insight_0.11.1      igraph_1.2.6        pkgconfig_2.0.3    
## [103] numDeriv_2016.8-1.1 foreign_0.8-81      binom_1.1-1        
## [106] xml2_1.3.2          dygraphs_1.1.1.6    estimability_1.3   
## [109] rvest_0.3.6         callr_3.5.1         digest_0.6.27      
## [112] parameters_0.10.1   fastGHQuad_1.0      rmarkdown_2.6      
## [115] cellranger_1.1.0    curl_4.3            shiny_1.5.0        
## [118] gtools_3.8.2        nloptr_1.2.2.2      lifecycle_0.2.0    
## [121] nlme_3.1-151        jsonlite_1.7.2      fansi_0.4.1        
## [124] pillar_1.4.7        loo_2.4.1           fastmap_1.0.1      
## [127] httr_1.4.2          plotrix_3.7-8       DEoptimR_1.0-8     
## [130] pkgbuild_1.2.0      survival_3.2-7      glue_1.4.2         
## [133] xts_0.12.1          bayestestR_0.8.0    zip_2.1.1          
## [136] shinythemes_1.1.2   iterators_1.0.13    stringi_1.5.3      
## [139] performance_0.6.1
# Read data ======================================
col_types <- cols(
  labid = col_factor(),
  subid = col_factor(),
  subid_unique = col_factor(),
  CDI.form = col_factor(),
  CDI.nwords = col_integer(),
  CDI.prop = col_number(),
  CDI.agerange = col_factor(),
  CDI.agedays = col_integer(),
  CDI.agemin = col_integer(),
  CDI.agemax = col_integer(),
  vocab_nwords = col_integer(),
  standardized.score.CDI = col_character(),
  standardized.score.CDI.num = col_number(),
  IDS_pref = col_number(),
  language = col_factor(),
  language_zone = col_factor(),
  CDI.error = col_logical(),
  Notes = col_character(),
  trial_order = col_factor(),
  method = col_factor(),
  age_days = col_integer(),
  age_mo = col_number(),
  age_group = col_factor(),
  nae = col_logical(),
  gender = col_factor(),
  second_session = col_logical()
)
data.total <- read_csv("data/02b_processed.csv", col_types = col_types)
# TODO: add saved results

Before moving on with the analysis, we have to ready the data by (a) checking for colinearity between z_age_months and CDI.z_age_months and correcting this if necessary, and (b) setting up the contrasts described in our data analysis.

1.1 Colinearity check

First, we run a Kappa test on the possibility of colinearity between z_age_months and CDI.z_age_months.

# Run kappa test
k.age_months <- model.matrix(~ z_age_months + CDI.z_age_months, data = data.total) %>%
  kappa(exact = T)

With a value of 3.1953332, we do not have a colinearity issue and can proceed with the data analysis as planned (The criteria of indicating colinearity is that kappa > 10).

1.2 Contrast Setups

We need gender as an effect-coded factor, and method as a deviation-coded factor. This is achieved in R by using the contr.sum() function with the number of levels for each factor. Notably, when subsetting the UK sample, only two levels of method out of the three in total were left.

# Set contrasts on the total dataset =============
contrasts(data.total$gender) <- contr.sum(2)
contrasts(data.total$method) <- contr.sum(3)

# Create sub-datasets, with contrasts ============
## NAE
data.nae <- data.total %>% subset(language_zone == "NAE") %>% droplevels()
contrasts(data.nae$gender) <- contr.sum(2)
contrasts(data.nae$method) <- contr.sum(3)

## UK
data.uk <- data.total %>% subset(language_zone == "British") %>% droplevels()
contrasts(data.uk$gender) <- contr.sum(2)
contrasts(data.uk$method) <- contr.sum(2) #note that UK sample has only 2 levels, so sum of zero contrasts set to 2 levels

## Other
data.other <- data.total %>% subset(language_zone == "Other") %>% droplevels()
contrasts(data.other$gender) <- contr.sum(2)
contrasts(data.other$method) <- contr.sum(3)

2 Descriptive Statistics

We first assess the amount of data we have overall per condition and their shape overall.

data.total %>%
  group_by(language_zone, CDI.agerange, method, gender) %>%
  summarise(N = n(), age = mean(CDI.agedays), sd = sd(CDI.agedays)) %>%
  kable()
## `summarise()` regrouping output by 'language_zone', 'CDI.agerange', 'method' (override with `.groups` argument)
language_zone CDI.agerange method gender N age sd
British 18 singlescreen F 24 551.2083 10.741950
British 18 singlescreen M 24 550.5000 13.497182
British 18 eyetracking F 9 548.5556 9.139353
British 18 eyetracking M 11 547.8182 10.147100
British 24 singlescreen F 18 741.6111 13.128948
British 24 singlescreen M 15 745.0667 15.549307
British 24 eyetracking F 13 731.5385 12.162890
British 24 eyetracking M 5 737.8000 8.228001
Other 18 singlescreen F 11 541.5455 7.160498
Other 18 singlescreen M 13 544.3077 6.663140
Other 18 eyetracking F 26 556.0769 14.133430
Other 18 eyetracking M 30 560.8333 16.128970
Other 18 hpp F 38 549.0526 10.812774
Other 18 hpp M 38 555.1579 13.664961
Other 24 singlescreen F 10 731.3000 14.802777
Other 24 singlescreen M 12 721.1667 13.888081
Other 24 eyetracking F 28 730.1786 9.996494
Other 24 eyetracking M 25 730.1600 7.711896
Other 24 hpp F 45 730.7333 17.357733
Other 24 hpp M 35 730.5143 15.849237
NAE 18 singlescreen F 19 554.9474 20.780530
NAE 18 singlescreen M 14 581.2143 24.925052
NAE 18 eyetracking F 1 549.0000 NA
NAE 18 hpp F 56 560.9821 21.584920
NAE 18 hpp M 57 559.6140 21.163272
NAE 24 singlescreen F 23 737.1739 26.604258
NAE 24 singlescreen M 20 739.6000 21.352678
NAE 24 eyetracking F 2 756.5000 34.648232
NAE 24 eyetracking M 1 745.0000 NA
NAE 24 hpp F 54 733.9630 27.476895
NAE 24 hpp M 60 727.9500 27.409899

More detailed information about Descriptive Statistics

#number of lab
data.total %>% 
  select(labid, language_zone) %>% 
  unique() %>% 
  group_by(language_zone) %>% 
  count()
## # A tibble: 3 x 2
## # Groups:   language_zone [3]
##   language_zone     n
##   <fct>         <int>
## 1 British           4
## 2 Other             8
## 3 NAE               8
data.total %>% 
  group_by(language_zone, CDI.agerange) %>% 
  summarize(N = n())
## `summarise()` regrouping output by 'language_zone' (override with `.groups` argument)
## # A tibble: 6 x 3
## # Groups:   language_zone [3]
##   language_zone CDI.agerange     N
##   <fct>         <fct>        <int>
## 1 British       18              68
## 2 British       24              51
## 3 Other         18             156
## 4 Other         24             155
## 5 NAE           18             147
## 6 NAE           24             160
# age range in each age group and language zone
data.total %>% 
  select(subid, language_zone, CDI.agedays, CDI.agerange) %>% 
  unique() %>% 
  group_by(language_zone, CDI.agerange) %>% 
  summarize(age_min = (min(CDI.agedays)/365.25*12),
            age_max = (max(CDI.agedays)/365.25*12))
## `summarise()` regrouping output by 'language_zone' (override with `.groups` argument)
## # A tibble: 6 x 4
## # Groups:   language_zone [3]
##   language_zone CDI.agerange age_min age_max
##   <fct>         <fct>          <dbl>   <dbl>
## 1 British       18              17.4    19.2
## 2 British       24              23.4    25.1
## 3 Other         18              17.4    19.5
## 4 Other         24              22.5    25.2
## 5 NAE           18              16.6    20.0
## 6 NAE           24              22.2    25.9

We then assess the data per lab in terms of sample size and CDI score (vocabulary size, for consistency between language zones).

by_lab <- data.total %>%
  group_by(labid, language_zone, CDI.agerange) %>%
  mutate(tested = n_distinct(subid_unique)) %>%
  select(labid, language_zone, CDI.agerange, tested, vocab_nwords) %>%
  nest(scores = vocab_nwords) %>%
  mutate(model = map(scores, ~ lm(vocab_nwords ~ 1, data = .x)),
         ci = map(model, confint)) %>%
  transmute(tested = tested,
            mean = map_dbl(model, ~ coefficients(.x)[[1]]),
            ci_lower = map_dbl(ci, 1),
            ci_upper = map_dbl(ci, 2)) %>%
  arrange(language_zone) %>%
  rownames_to_column()

# TODO: find a way to group by language zone?
ggplot(by_lab,
       aes(x = labid, colour = language_zone,
           y = mean, ymin = ci_lower, ymax = ci_upper)) + 
  geom_linerange() + 
  geom_point(aes(size = tested)) + 
  facet_grid(cols = vars(CDI.agerange), scales = "free") + coord_flip(ylim = c(0, 500)) +
  xlab("Lab") + ylab("Vocabulary size") +
  scale_colour_brewer(palette = "Dark2", name = "Language zone",
                      breaks = c("British", "NAE", "Other")) +
  scale_size_continuous(name = "N", breaks = function(x) c(min(x), mean(x), max(x))) +
  theme(legend.position = "bottom")

3 Sample Theory Based Statistics

3.1 Simple Correlation

First, we want to assess quickly if there is a direct correlation between IDS preference and CDI score, computing a Pearson#’s product-moment correlation. We use standardized CDI scores for the North American sample, and raw scores for the British sample. [NOTE FROM MELANIE: I’m not sure this makes sense to do with the raw UK scores - since we’re collapsing across 18 and 24 month data and the data aren#’t standardized by age.].

# Statistics =====================================
## North American Sample
test.pearson.nae <- cor.test(data.nae$IDS_pref,
                             data.nae$z_standardized_CDI,
                             alternative = "two.sided", method = "pearson")

test.pearson.nae
## 
##  Pearson's product-moment correlation
## 
## data:  data.nae$IDS_pref and data.nae$z_standardized_CDI
## t = -0.91847, df = 305, p-value = 0.3591
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.16349833  0.05977293
## sample estimates:
##         cor 
## -0.05251901
## UK Sample
test.pearson.uk <- cor.test(data.uk$IDS_pref,
                            data.uk$z_vocab_nwords,
                            alternative = "two.sided", method = "pearson")

test.pearson.uk
## 
##  Pearson's product-moment correlation
## 
## data:  data.uk$IDS_pref and data.uk$z_vocab_nwords
## t = 1.0327, df = 117, p-value = 0.3039
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08643398  0.27040987
## sample estimates:
##        cor 
## 0.09504018

Plots for correlation

## North American Sample
### Get correlation value for annotation
cor_text <- "paste(italic(R)^2, \" =\")"
cor_value <- round(test.pearson.nae$estimate, 3)

### Build plot
plot.pearson.nae <- data.nae %>%
  ggplot(aes(x = IDS_pref,
             y = standardized.score.CDI.num)) +
  xlab("IDS preference") + ylab("Standardized CDI score") +
  geom_point() +
  geom_smooth(method = lm) +
  annotate("text", x = -.9, y = 51, parse = T, size = 4,
           label = paste(cor_text, cor_value, sep = "~"))

## UK Sample
cor_value <- round(test.pearson.uk$estimate, 3)
plot.pearson.uk <- data.uk %>%
  ggplot(aes(x = IDS_pref,
             y = vocab_nwords)) +
  xlab("IDS preference") + ylab("Vocabulary size (in words)") +
  geom_point() +
  geom_smooth(method = lm) +
  annotate("text", x = .8, y = 150, parse = T, size = 4,
           label = paste(cor_text, cor_value, sep = "~"))

# Global plot
plot.pearson <- ggarrange(plot.pearson.nae, plot.pearson.uk, ncol = 2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

plot.pearson

ggsave("plots/corr_plot.pdf", plot.pearson,
       units = "mm", width = 180, height = 100, dpi = 1000)

We see no obvious direct link between IDS prefernce and CDI score here. However, an effect might appear once we take into account various factors that might interact with IDS preference and/or CDI score. We can also first enhance these plots with information about the age group at which infants were tested (18- or 24-month-old), using vocabulary size to better compare the NAE and UK samples.

plot.age_group <- data.total %>%
  subset(language_zone != "Other") %>%
  droplevels() %>%
  ggplot(aes(x = IDS_pref,
             y = vocab_nwords,
             colour = CDI.agerange)) +
  facet_wrap(vars(language_zone),
             labeller = as_labeller(c("British" = "UK samples",
                                      "NAE" = "North Amercian samples"))) +
  xlab("Standardized IDS prefence score") + ylab("Vocabulary size (in words)") + theme(legend.position = "top") +
  geom_point() +
  geom_smooth(method = lm) +
  scale_colour_brewer(palette = "Dark2", name = "Age group",
                      breaks = c("18", "24"),
                      labels = c("18mo", "24m"))
ggsave("plots/scatter_age.pdf", plot.age_group,
       units = "mm", width = 180, height = 100, dpi = 1000)
## `geom_smooth()` using formula 'y ~ x'
(plot.age_group)
## `geom_smooth()` using formula 'y ~ x'

3.2 Mixed-Effects Model by Language Zone

Here, we run a mixed-effects model including only theoretically motivated effects, as described in the pre-registration. We start with the full model bellow, simplifying the random effects structure until it converges.

3.2.1 NAE full model

# Run models =====================================
## NAE

data.nae$centered_IDS_pref <- scale(data.nae$IDS_pref, center = T, scale = F)

lmer.full.nae <- lmer(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + centered_IDS_pref +
                        centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months +
                        (1 | labid),
                      data = data.nae)

summary(lmer.full.nae)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months +  
##     method + centered_IDS_pref + centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months +  
##     centered_IDS_pref:z_age_months + (1 | labid)
##    Data: data.nae
## 
## REML criterion at convergence: 2861
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7649 -0.8396 -0.1049  0.7743  2.3582 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  labid    (Intercept)  78.28    8.847  
##  Residual             752.79   27.437  
## Number of obs: 307, groups:  labid, 8
## 
## Fixed effects:
##                                     Estimate Std. Error        df t value
## (Intercept)                         31.38428    6.13351  28.72420   5.117
## CDI.z_age_months                     0.07662    0.55332 295.20197   0.138
## gender1                             -1.07741    1.61552 291.53707  -0.667
## z_age_months                         0.35211    0.56449 153.40385   0.624
## method1                              8.62370    5.88516  96.93942   1.465
## method2                            -16.48345   10.24023 225.60947  -1.610
## centered_IDS_pref                  -36.37414   23.08356 289.74626  -1.576
## method1:centered_IDS_pref           27.84822   23.61831 290.27716   1.179
## method2:centered_IDS_pref          -58.54287   46.03821 289.28918  -1.272
## CDI.z_age_months:centered_IDS_pref   1.15753    1.67303 290.78963   0.692
## z_age_months:centered_IDS_pref       1.78245    1.73910 292.19209   1.025
##                                    Pr(>|t|)    
## (Intercept)                        1.88e-05 ***
## CDI.z_age_months                      0.890    
## gender1                               0.505    
## z_age_months                          0.534    
## method1                               0.146    
## method2                               0.109    
## centered_IDS_pref                     0.116    
## method1:centered_IDS_pref             0.239    
## method2:centered_IDS_pref             0.205    
## CDI.z_age_months:centered_IDS_pref    0.490    
## z_age_months:centered_IDS_pref        0.306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) CDI.z__ gendr1 z_g_mn methd1 methd2 c_IDS_ m1:_ID m2:_ID
## CDI.z_g_mnt -0.079                                                         
## gender1     -0.063  0.031                                                  
## z_age_mnths  0.031 -0.043  -0.070                                          
## method1     -0.562 -0.003   0.019  0.140                                   
## method2      0.800 -0.045  -0.057 -0.010 -0.777                            
## cntrd_IDS_p  0.202  0.019   0.027 -0.037 -0.203  0.238                     
## mthd1:_IDS_ -0.184 -0.049  -0.014  0.019  0.231 -0.244 -0.929              
## mthd2:_IDS_  0.195  0.040   0.047 -0.012 -0.222  0.244  0.970 -0.972       
## CDI.__:_IDS -0.012 -0.021  -0.030 -0.018 -0.043  0.018 -0.075  0.038 -0.072
## z_g_m:_IDS_ -0.029  0.044   0.131  0.035 -0.049  0.030  0.037 -0.063  0.121
##             CDI.__:
## CDI.z_g_mnt        
## gender1            
## z_age_mnths        
## method1            
## method2            
## cntrd_IDS_p        
## mthd1:_IDS_        
## mthd2:_IDS_        
## CDI.__:_IDS        
## z_g_m:_IDS_ -0.097
robust_lmer.full.nae <- robustlmm::rlmer(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + centered_IDS_pref +
                        centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months +
                        (1 | labid),
                      data = data.nae)


summary(lmer.full.nae)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months +  
##     method + centered_IDS_pref + centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months +  
##     centered_IDS_pref:z_age_months + (1 | labid)
##    Data: data.nae
## 
## REML criterion at convergence: 2861
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7649 -0.8396 -0.1049  0.7743  2.3582 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  labid    (Intercept)  78.28    8.847  
##  Residual             752.79   27.437  
## Number of obs: 307, groups:  labid, 8
## 
## Fixed effects:
##                                     Estimate Std. Error        df t value
## (Intercept)                         31.38428    6.13351  28.72420   5.117
## CDI.z_age_months                     0.07662    0.55332 295.20197   0.138
## gender1                             -1.07741    1.61552 291.53707  -0.667
## z_age_months                         0.35211    0.56449 153.40385   0.624
## method1                              8.62370    5.88516  96.93942   1.465
## method2                            -16.48345   10.24023 225.60947  -1.610
## centered_IDS_pref                  -36.37414   23.08356 289.74626  -1.576
## method1:centered_IDS_pref           27.84822   23.61831 290.27716   1.179
## method2:centered_IDS_pref          -58.54287   46.03821 289.28918  -1.272
## CDI.z_age_months:centered_IDS_pref   1.15753    1.67303 290.78963   0.692
## z_age_months:centered_IDS_pref       1.78245    1.73910 292.19209   1.025
##                                    Pr(>|t|)    
## (Intercept)                        1.88e-05 ***
## CDI.z_age_months                      0.890    
## gender1                               0.505    
## z_age_months                          0.534    
## method1                               0.146    
## method2                               0.109    
## centered_IDS_pref                     0.116    
## method1:centered_IDS_pref             0.239    
## method2:centered_IDS_pref             0.205    
## CDI.z_age_months:centered_IDS_pref    0.490    
## z_age_months:centered_IDS_pref        0.306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) CDI.z__ gendr1 z_g_mn methd1 methd2 c_IDS_ m1:_ID m2:_ID
## CDI.z_g_mnt -0.079                                                         
## gender1     -0.063  0.031                                                  
## z_age_mnths  0.031 -0.043  -0.070                                          
## method1     -0.562 -0.003   0.019  0.140                                   
## method2      0.800 -0.045  -0.057 -0.010 -0.777                            
## cntrd_IDS_p  0.202  0.019   0.027 -0.037 -0.203  0.238                     
## mthd1:_IDS_ -0.184 -0.049  -0.014  0.019  0.231 -0.244 -0.929              
## mthd2:_IDS_  0.195  0.040   0.047 -0.012 -0.222  0.244  0.970 -0.972       
## CDI.__:_IDS -0.012 -0.021  -0.030 -0.018 -0.043  0.018 -0.075  0.038 -0.072
## z_g_m:_IDS_ -0.029  0.044   0.131  0.035 -0.049  0.030  0.037 -0.063  0.121
##             CDI.__:
## CDI.z_g_mnt        
## gender1            
## z_age_mnths        
## method1            
## method2            
## cntrd_IDS_p        
## mthd1:_IDS_        
## mthd2:_IDS_        
## CDI.__:_IDS        
## z_g_m:_IDS_ -0.097
full.nae_pvalue <- anova(lmer.full.nae) %>% 
  as_tibble(rownames = "Parameter") #this gives us the Type III p values

# ==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref
# method
# z_age_months
# gender
# CDI.z_age_months
#==========

3.2.1.1 (Optional) Checking mixed-model assumptions. We will check the following:

  1. Linearlity
  2. Normality of the residuals
  3. Homoscedasticity of residuals
  4. No autocorrelation
  5. No multicollinearity
#First, check linearity
# data.nae$resid <- residuals(lmer.full.nae)
# 
# plot(data.nae$resid, data.nae$standardized.score.CDI)

#Second, check normality
plot_model(lmer.full.nae, type = 'diag') ## we do have right-skewed normality of residuals
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.18
## Current Matrix version is 1.3.2
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

#Third, check autocorrelation
re_run_lme.full.nae <- nlme::lme(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + centered_IDS_pref +
                        centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months, random = ~1 | labid,
                        method = "REML",
                      data = data.nae, na.action = na.exclude)

plot(nlme::ACF(re_run_lme.full.nae, resType = "normalized")) #there is no sign for autocorrelation

#Lastly, check multi-collinearity
car::vif(lmer.full.nae) #we do see a multicollineartiy for the IDS preference variable, even though we have centered the IDS preference score. It is probably related to the number the participating labs (as this is the group level that we are controlling) and how we entered interaction between IDS preference and other variables (that lack variability in the current sample). We need to keep IDS preference in the model as exploring the relationship between IDS preference and CDI score is the key research question in the paper.
##                                         GVIF Df GVIF^(1/(2*Df))
## CDI.z_age_months                    1.018050  1        1.008985
## gender                              1.052996  1        1.026156
## z_age_months                        1.076038  1        1.037322
## method                              1.147009  2        1.034884
## centered_IDS_pref                  22.501425  1        4.743567
## method:centered_IDS_pref           25.126190  2        2.238884
## CDI.z_age_months:centered_IDS_pref  1.032456  1        1.016099
## z_age_months:centered_IDS_pref      1.306986  1        1.143235

3.2.2 UK full model

lmer.full.uk <- lmer(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                       (1 | labid),
                     data = data.uk) 

summary(lmer.full.uk)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method +  
##     IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months +  
##     IDS_pref:z_age_months + (1 | labid)
##    Data: data.uk
## 
## REML criterion at convergence: 1342.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0447 -0.6691 -0.0987  0.5147  3.5513 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  labid    (Intercept)  148.7   12.20   
##  Residual             7856.7   88.64   
## Number of obs: 119, groups:  labid, 4
## 
## Fixed effects:
##                           Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                175.103     11.540   3.162  15.174 0.000464 ***
## CDI.z_age_months            28.628      2.824 109.795  10.138  < 2e-16 ***
## gender1                     17.578      8.623 109.015   2.039 0.043906 *  
## z_age_months                -6.034      3.042 101.094  -1.984 0.050013 .  
## method1                    -16.483     12.871   4.413  -1.281 0.263499    
## IDS_pref                     4.201     23.986 109.837   0.175 0.861288    
## method1:IDS_pref             4.413     27.446 109.729   0.161 0.872550    
## CDI.z_age_months:IDS_pref    0.600      7.877 108.111   0.076 0.939422    
## z_age_months:IDS_pref        1.143      8.291 109.890   0.138 0.890629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) CDI.z__ gendr1 z_g_mn methd1 IDS_pr m1:IDS CDI.__:
## CDI.z_g_mnt  0.148                                                   
## gender1     -0.023 -0.099                                            
## z_age_mnths -0.005 -0.106  -0.001                                    
## method1     -0.287 -0.075  -0.029  0.454                             
## IDS_pref    -0.314 -0.092  -0.118  0.059  0.179                      
## mthd1:IDS_p  0.163  0.122  -0.012 -0.119 -0.264 -0.157               
## CDI.__:IDS_ -0.126 -0.328   0.076  0.063  0.103  0.112 -0.195        
## z_g_mn:IDS_  0.014  0.095  -0.232 -0.360 -0.106 -0.096  0.493 -0.280
full.uk_pvalue <- anova(lmer.full.uk) %>% 
  as_tibble(rownames = "Parameter") #this gives us the Type III p values

#==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref
# method
# z_age_months
# gender
# CDI.z_age_months

3.2.2.1 (Optional) Checking mixed-model assumptions. We will check the following:

  1. Linearlity
  2. Normality of the residuals
  3. Homoscedasticity of residuals
  4. No autocorrelation
  5. No multicollinearity
#First, check linearity. The plot looks linear
data.uk$resid <- residuals(lmer.full.uk)
 
plot(data.uk$resid, data.uk$vocab_nwords)

#Second, check normality
plot_model(lmer.full.uk, type = 'diag') ## we do have right-skewed normality of residuals
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

#Third, check autocorrelation
re_run_lme.full.uk <- nlme::lme(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months, random = ~1 | labid,
                        method = "REML",
                      data = data.nae, na.action = na.exclude)

plot(nlme::ACF(re_run_lme.full.uk, resType = "normalized")) #there is no sign for autocorrelation

#Lastly, check multi-collinearity
car::vif(lmer.full.uk) #no problem for multicollinearlity
##          CDI.z_age_months                    gender              z_age_months 
##                  1.155355                  1.118400                  1.553116 
##                    method                  IDS_pref           method:IDS_pref 
##                  1.445581                  1.078478                  1.515016 
## CDI.z_age_months:IDS_pref     z_age_months:IDS_pref 
##                  1.227720                  1.819286

We now want to check the statistical power of significant effects, and discard any models with significant effects that do not reach 80% power. This however leads to too many warnings of singularity issues on the model updates inherent to the simr power simulations, hence we cannot obtain satisfactory power estimates as pre-registered.

AST: Note that we don’t have any IV(s) that turned out to be significant in the Full NAE model. So we won’t run the power analysis check. For the UK full model, there are two statistically significant IV: CDI_age and gender. The post hoc power check suggested that we have high power in detecting the effect of CDI_age but not gender. Note that gender has a smaller effect size to begin with, so this may partially explain why we have less power in detecting it in the model. As there can be a number of different factors that determines the posthoc power, we decided not to remove gender in the model based on posthoc power analysis check.

check_pwr_uk_cdi_age <- simr::powerSim(lmer.full.uk, test = fixed("CDI.z_age_months", method = "z"), seed = 2, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_uk_cdi_age

check_pwr_uk_gender <- simr::powerSim(lmer.full.uk, test = fixed("gender", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_uk_gender

3.2.3 Combined Sample

For this combined analysis, we first need to restrain the age range for the NAE sample (previously ±2 months, now ±1.5 months).

# Create dataset with British and NAE only
data.uk_nae <- data.total %>%
  subset(language_zone %in% c("British", "NAE")) %>%
  mutate(CDI.agemin = ifelse(language_zone == "NAE",
                             CDI.agemin + round(.5*365.25/12),
                             CDI.agemin),
         CDI.agemax = ifelse(language_zone == "NAE",
                             CDI.agemax - round(.5*365.25/12),
                             CDI.agemax)) %>%
  subset(!(CDI.agedays < CDI.agemin | CDI.agedays > CDI.agemax)) %>%
  droplevels()
# Create contrasts for analysis
contrasts(data.uk_nae$gender) <- contr.sum(2)
contrasts(data.uk_nae$method) <- contr.sum(3)
contrasts(data.uk_nae$language_zone) <- contr.sum(2)

We can then run the planned combined analysis adding the main effect and interactions of language_zone.

lmer.full.uk_nae <- lmer(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + IDS_pref:language_zone +
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                           (1 + CDI.z_age_months | labid),
                         data = data.uk_nae)

summary(lmer.full.uk_nae)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months +  
##     method + IDS_pref + IDS_pref:language_zone + IDS_pref:method +  
##     IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 +  
##     CDI.z_age_months | labid)
##    Data: data.uk_nae
## 
## REML criterion at convergence: -16.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2182 -0.6666 -0.1781  0.5001  3.4898 
## 
## Random effects:
##  Groups   Name             Variance  Std.Dev. Corr
##  labid    (Intercept)      0.0027125 0.05208      
##           CDI.z_age_months 0.0008073 0.02841  0.40
##  Residual                  0.0439990 0.20976      
## Number of obs: 401, groups:  labid, 12
## 
## Fixed effects:
##                             Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 0.323365   0.022302   6.812745  14.499 2.26e-06 ***
## CDI.z_age_months            0.042960   0.009300  12.434234   4.619 0.000539 ***
## language_zone1              0.080976   0.027480   7.098782   2.947 0.021166 *  
## gender1                     0.037004   0.010822 372.652269   3.419 0.000697 ***
## z_age_months               -0.001733   0.003994 289.264152  -0.434 0.664725    
## method1                    -0.010715   0.027902  12.016228  -0.384 0.707678    
## method2                    -0.001692   0.042360  10.268417  -0.040 0.968902    
## IDS_pref                   -0.005438   0.038612 382.246227  -0.141 0.888062    
## language_zone1:IDS_pref     0.068171   0.053566 374.217018   1.273 0.203927    
## method1:IDS_pref           -0.018179   0.051697 385.237670  -0.352 0.725301    
## method2:IDS_pref           -0.049967   0.080989 382.643250  -0.617 0.537625    
## CDI.z_age_months:IDS_pref   0.008436   0.011007 383.996827   0.766 0.443883    
## z_age_months:IDS_pref      -0.003172   0.011006 376.867698  -0.288 0.773312    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
combined.full.uk_nae_pvalue <- anova(lmer.full.uk_nae) %>% 
  as_tibble(rownames = "Parameter") #this gives us the Type III p values

#==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref:language_zone
# IDS_pref
# method
# z_age_months
# gender
# language_zone
#==========

3.2.3.1 (Optional) Checking mixed-model assumptions

  1. Linearlity
  2. Normality of the residuals
  3. Homoscedasticity of residuals
  4. No autocorrelation
  5. No multicollinearity
#First, check linearity. The plot looks linear
data.uk_nae$resid <- residuals(lmer.full.uk_nae)
 
plot(data.uk_nae$resid, data.uk_nae$CDI.prop)

#Second, check normality
plot_model(lmer.full.uk_nae, type = 'diag') ## we do have right-skewed normality of residuals
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

#Third, check autocorrelation
re_run_lme.full.uk_nae <- nlme::lme(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + IDS_pref:language_zone +
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months,
                           random = ~CDI.z_age_months  | labid,
                        method = "REML",
                      data = data.uk_nae, na.action = na.exclude)

plot(nlme::ACF(re_run_lme.full.uk_nae, resType = "normalized")) #there is no sign for autocorrelation

#Lastly, check multi-collinearity
car::vif(lmer.full.uk_nae) #no problem for multicollinearlity
##                               GVIF Df GVIF^(1/(2*Df))
## CDI.z_age_months          1.032963  1        1.016348
## language_zone             1.979191  1        1.406837
## gender                    1.039619  1        1.019617
## z_age_months              1.328619  1        1.152657
## method                    2.447762  2        1.250813
## IDS_pref                  1.449858  1        1.204101
## language_zone:IDS_pref    2.900147  1        1.702982
## method:IDS_pref           3.854846  2        1.401205
## CDI.z_age_months:IDS_pref 1.078314  1        1.038419
## z_age_months:IDS_pref     1.509419  1        1.228584

We then compute \(p\)-values, but leave out power estimates for those \(p\)-values as above. Again, we have a lot of singular fit issues for the power checks and decided not to remove parameters based on posthoc power analysis.

check_pwr_combined_cdi_age <- simr::powerSim(lmer.full.uk_nae, test = fixed("CDI.z_age_months", method = "z"), seed = 2, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_combined_cdi_age

check_pwr_combined_lang_zone <- simr::powerSim(lmer.full.uk_nae, test = fixed("language_zone", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_combined_lang_zone

check_pwr_combined_gender <- simr::powerSim(lmer.full.uk_nae, test = fixed("gender", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_combined_gender

4 Bayesian Statistics

Given that the `frequentist#’ tests above did not discover any expected significant effects, we now run a Bayesian analysis, specifically Bayes factors (or \(b\)-values) from model comparisons, to determine whether our data provide evidence for a true null effect or merely fail to provide evidence either way.

We first run models on the separate UK and NAE samples.

First of all, we need to set up the priors for the Bayeisan analysis. In the RR, we said that we would use a Cauchy distribution for priors over fixed effects with normally distributed priors over random effect parameters."

4.1 Set up priors

Set up a really weak and not too informative prior

prior <-c(set_prior("normal(10, 1)", class = "Intercept"), #we can't have negative CDI
              set_prior("cauchy(0, 1)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd")) 

prior_combined_mod <-c(set_prior("normal(10, 1)", class = "Intercept"),
              set_prior("cauchy(0, 1)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

Set up another prior that mainly changed the meaan for the fixed effect and see if there is big difference? Note: numbers changed but not to the extent that change BF conclusion.

prior_big <-c(set_prior("normal(10, 1)", class = "Intercept"), #we can't have negative CDI
              set_prior("cauchy(10, 5)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd")) 

prior_combined_big <-c(set_prior("normal(10, 1)", class = "Intercept"),
              set_prior("cauchy(10, 5)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

We also try more informative prior.

4.2 NAE info prior

prior_nae_1 <-c(set_prior("normal(50, 1)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               set_prior("normal(0, 1)", class = "sd")) 

prior_nae_null_1 <-c(set_prior("normal(50, 1)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               set_prior("normal(0, 1)", class = "sd")) 

prior_nae_2 <-c(set_prior("normal(50, 1)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
              #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_nae_null_2 <-c(set_prior("normal(50, 1)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

4.3 UK info prior

prior_uk_1 <-c(set_prior("normal(40, 1)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_uk_null_1 <-c(set_prior("normal(40, 1)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_uk_2 <-c(set_prior("normal(40, 1)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_uk_null_2 <-c(set_prior("normal(40, 1)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

4.4 Combined info prior

# Note that NAE CDI WS total score is 680 and UK CDI total score is 416. 

# NAE average scores for 50th percentile is 76.3, so it is 76.3/680 ~ 11% and UK CDI average scores for 50th percentile is 41.6, so it is 41.6/416 ~ 10%

prior_combined_full_info <-c(set_prior("normal(10, 1)", class = "Intercept"), # given the above info, we expect on average, proportion score is 10% 
              set_prior("cauchy(3, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase in UK cdi and lead to 28 words increase in NAE (proxy from wordbank). In general, 15/416 = 3% and 28/680 = 4% 
              set_prior("cauchy(-3, 1)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males, 11/416 ~ 2.6% and Wordbank NAE suggested females knows approximately 27 words than males, 27/680 ~ 3.9%
               set_prior("cauchy(0, 1)", class = "b", coef = "language_zone1"), #weak prior as don't know the effect of language zone to CDI score 
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "language_zone1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

prior_combined_null_info <-c(set_prior("normal(10, 1)", class = "Intercept"), # given the above info, we expect on average, proportion score is 10% 
              set_prior("cauchy(3, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase in UK cdi and lead to 28 words increase in NAE (proxy from wordbank). In general, 15/416 = 3% and 28/680 = 4% 
              set_prior("cauchy(-3, 1)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males, 11/416 ~ 2.6% and Wordbank NAE suggested females knows approximately 27 words than males, 27/680 ~ 3.9%
               set_prior("cauchy(0, 1)", class = "b", coef = "language_zone1"), #weak prior as don't know the effect of language zone to CDI score 
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

In the following analysis, we follow what we have proposed in the RR about calculating the Bayes factor for the three variables of interest: 1) The main effect IDS preference on CDI scores 2) The interaction effect between IDS preference and IDS test age on CDI scores 3) The interaction effect between language_zone (i.e., dialect in RR) and IDS preference on CDI scores

For the first two variables, we run it on the NAE and UK model separately. For the third variable, we run in on the combined UK & NAE model.

To calcualte Bayes factor for the main effect IDS preference on CDI scores, we need to two models: a full model that contains the main effect of variables of interest, and the other null model that does not contain the main effect of interest. As all the interaction terms in the models contains IDS preference, in the following, we only include main effects in the full and null models for comparisons.

4.5 NAE model (main effect of IDS preference)

Sys.setenv(R_FUTURE_RNG_ONMISUSE = "ignore") #this line of code is needed to get rid of warning message due to the "future" package. https://github.com/satijalab/seurat/issues/3622 

lmer.bf_full_1.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + IDS_pref +
                        (1 | labid),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=456)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '3b2c8dc852b93f8f4a16392e0b6193a9' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.742789 seconds (Warm-up)
## Chain 1:                0.769977 seconds (Sampling)
## Chain 1:                1.51277 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '3b2c8dc852b93f8f4a16392e0b6193a9' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.741755 seconds (Warm-up)
## Chain 2:                0.472661 seconds (Sampling)
## Chain 2:                1.21442 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '3b2c8dc852b93f8f4a16392e0b6193a9' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 5.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.54 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.69613 seconds (Warm-up)
## Chain 3:                0.749866 seconds (Sampling)
## Chain 3:                1.446 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '3b2c8dc852b93f8f4a16392e0b6193a9' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 4.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.693584 seconds (Warm-up)
## Chain 4:                0.662766 seconds (Sampling)
## Chain 4:                1.35635 seconds (Total)
## Chain 4:
summary(lmer.bf_full_1.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + (1 | labid) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.18      1.07     0.17     4.19 1.00     2834     2875
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           47.60      1.67    44.19    50.82 1.00     6735     4199
## CDI.z_age_months     0.15      0.57    -0.99     1.28 1.00     9418     5836
## gender1             -1.98      1.69    -5.38     1.24 1.00     9021     6193
## z_age_months        -0.04      0.44    -0.88     0.82 1.00     8768     5974
## method1             -0.01      1.30    -2.70     2.72 1.00     6826     4206
## method2             -0.28      1.97    -4.88     3.47 1.00     6414     3069
## IDS_pref            -0.57      2.04    -5.98     2.79 1.00     7243     2554
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.66      1.25    26.34    31.23 1.00     8001     5743
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_1.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + 
                        (1 | labid),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_null_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =123)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '32a715b1500064463a75f1c5a10e4ab8' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5.9e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.657162 seconds (Warm-up)
## Chain 1:                0.438973 seconds (Sampling)
## Chain 1:                1.09614 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '32a715b1500064463a75f1c5a10e4ab8' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.9 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.71534 seconds (Warm-up)
## Chain 2:                0.801815 seconds (Sampling)
## Chain 2:                1.51716 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '32a715b1500064463a75f1c5a10e4ab8' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 8.5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.85 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.58087 seconds (Warm-up)
## Chain 3:                0.812227 seconds (Sampling)
## Chain 3:                1.3931 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '32a715b1500064463a75f1c5a10e4ab8' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 9.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.93 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.706565 seconds (Warm-up)
## Chain 4:                0.890216 seconds (Sampling)
## Chain 4:                1.59678 seconds (Total)
## Chain 4:
summary(lmer.bf_null_1.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + (1 | labid) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.18      1.07     0.17     4.14 1.00     3174     2530
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           47.50      1.63    44.04    50.60 1.00     8004     4014
## CDI.z_age_months     0.13      0.57    -0.98     1.24 1.00    11260     5283
## gender1             -1.96      1.67    -5.24     1.30 1.00    12196     5268
## z_age_months        -0.05      0.43    -0.90     0.80 1.00    11669     5684
## method1              0.00      1.30    -2.72     2.69 1.00     8248     4403
## method2             -0.30      1.94    -4.93     3.34 1.00     6939     3312
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.66      1.25    26.33    31.19 1.00     8247     5849
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_1.nae <- bridge_sampler(lmer.bf_full_1.nae, silent = T)
marloglik_null_1.nae <- bridge_sampler(lmer.bf_null_1.nae, silent = T)

Bayes factor

BF_full_1.nae <- bayes_factor(marloglik_full_1.nae, marloglik_null_1.nae)
BF_full_1.nae
## Estimated Bayes factor in favor of x1 over x2: 0.91418

4.6 NAE model (interaction between IDS test age and IDS preference)

lmer.bf_full_2.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + IDS_pref +
                        IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                        (1 | labid),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=890)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'c92e15777c97d2cdb32465f3d1c85d38' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.74 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.987588 seconds (Warm-up)
## Chain 1:                0.771481 seconds (Sampling)
## Chain 1:                1.75907 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'c92e15777c97d2cdb32465f3d1c85d38' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7.8e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.802676 seconds (Warm-up)
## Chain 2:                0.863793 seconds (Sampling)
## Chain 2:                1.66647 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'c92e15777c97d2cdb32465f3d1c85d38' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6.9e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.69 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.812938 seconds (Warm-up)
## Chain 3:                0.726661 seconds (Sampling)
## Chain 3:                1.5396 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'c92e15777c97d2cdb32465f3d1c85d38' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 7.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.72 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.868436 seconds (Warm-up)
## Chain 4:                0.790037 seconds (Sampling)
## Chain 4:                1.65847 seconds (Total)
## Chain 4:
summary(lmer.bf_full_2.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 | labid) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.15      1.05     0.17     4.06 1.00     3214     3025
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    47.44      1.74    43.88    50.66 1.00     7073
## CDI.z_age_months              0.08      0.60    -1.08     1.24 1.00     8998
## gender1                      -1.90      1.66    -5.10     1.34 1.00     9048
## z_age_months                 -0.12      0.45    -1.00     0.77 1.00     8524
## method1                      -0.01      1.35    -2.81     2.73 1.00     7236
## method2                      -0.34      2.06    -5.14     3.45 1.00     6897
## IDS_pref                     -0.84      2.38    -7.36     2.66 1.00     6695
## method1:IDS_pref             -0.00      2.00    -4.42     4.30 1.00     8656
## method2:IDS_pref              0.28      2.36    -4.27     6.11 1.00     7900
## CDI.z_age_months:IDS_pref     0.44      1.10    -1.55     2.94 1.00     8479
## z_age_months:IDS_pref         0.76      1.15    -1.16     3.45 1.00     7334
##                           Tail_ESS
## Intercept                     4191
## CDI.z_age_months              5991
## gender1                       5386
## z_age_months                  6523
## method1                       3651
## method2                       3111
## IDS_pref                      2868
## method1:IDS_pref              3268
## method2:IDS_pref              2693
## CDI.z_age_months:IDS_pref     4637
## z_age_months:IDS_pref         4296
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.66      1.26    26.32    31.19 1.00     7384     6053
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_2.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + IDS_pref +
                        IDS_pref:method + IDS_pref:CDI.z_age_months + 
                        (1 | labid),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_null_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =102)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '91b614181d9419722a0fc80a47971299' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000262 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.62 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.83679 seconds (Warm-up)
## Chain 1:                0.982422 seconds (Sampling)
## Chain 1:                1.81921 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '91b614181d9419722a0fc80a47971299' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.74 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.871272 seconds (Warm-up)
## Chain 2:                0.675621 seconds (Sampling)
## Chain 2:                1.54689 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '91b614181d9419722a0fc80a47971299' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.962901 seconds (Warm-up)
## Chain 3:                0.55674 seconds (Sampling)
## Chain 3:                1.51964 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '91b614181d9419722a0fc80a47971299' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 7.6e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.76 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.905299 seconds (Warm-up)
## Chain 4:                0.732245 seconds (Sampling)
## Chain 4:                1.63754 seconds (Total)
## Chain 4:
summary(lmer.bf_null_2.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + (1 | labid) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.19      1.06     0.18     4.15 1.00     2544     2153
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    47.54      1.71    43.86    50.74 1.00     6544
## CDI.z_age_months              0.07      0.60    -1.09     1.25 1.00     7444
## gender1                      -1.99      1.68    -5.33     1.30 1.00     9262
## z_age_months                 -0.04      0.43    -0.91     0.81 1.00     7060
## method1                       0.01      1.32    -2.85     2.90 1.00     6039
## method2                      -0.37      1.98    -5.33     3.15 1.00     6133
## IDS_pref                     -0.66      2.29    -6.80     3.04 1.00     5705
## method1:IDS_pref             -0.13      1.99    -4.73     4.07 1.00     7113
## method2:IDS_pref              0.13      2.39    -4.60     5.66 1.00     5917
## CDI.z_age_months:IDS_pref     0.47      1.13    -1.56     3.03 1.00     7947
##                           Tail_ESS
## Intercept                     3678
## CDI.z_age_months              5504
## gender1                       5942
## z_age_months                  5475
## method1                       3608
## method2                       2658
## IDS_pref                      2326
## method1:IDS_pref              3277
## method2:IDS_pref              2471
## CDI.z_age_months:IDS_pref     4550
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.65      1.23    26.40    31.20 1.00     7405     6208
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_2.nae <- bridge_sampler(lmer.bf_full_2.nae, silent = T)
marloglik_null_2.nae <- bridge_sampler(lmer.bf_null_2.nae, silent = T)

Bayes factor

BF_full_2.nae <- bayes_factor(marloglik_full_2.nae, marloglik_null_2.nae)
BF_full_2.nae
## Estimated Bayes factor in favor of x1 over x2: 0.81012

4.7 UK model (main effect of IDS preference)

lmer.bf_full_1.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                      (1 | labid),
                     data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=213)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '06bfb59652a0c5924487daf6d3cbc793' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 6.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.622256 seconds (Warm-up)
## Chain 1:                0.476119 seconds (Sampling)
## Chain 1:                1.09838 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '06bfb59652a0c5924487daf6d3cbc793' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.78034 seconds (Warm-up)
## Chain 2:                0.514808 seconds (Sampling)
## Chain 2:                1.29515 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '06bfb59652a0c5924487daf6d3cbc793' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7.8e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.679524 seconds (Warm-up)
## Chain 3:                0.907757 seconds (Sampling)
## Chain 3:                1.58728 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '06bfb59652a0c5924487daf6d3cbc793' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000222 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.22 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.740539 seconds (Warm-up)
## Chain 4:                0.838729 seconds (Sampling)
## Chain 4:                1.57927 seconds (Total)
## Chain 4:
summary(lmer.bf_full_1.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + (1 | labid) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.84      0.63     0.04     2.35 1.00     5216     2823
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           49.61      3.35    43.36    56.17 1.00     3288     3757
## CDI.z_age_months    22.88      5.82    14.39    34.27 1.00     3814     5617
## gender1              0.83     12.22   -15.11    30.57 1.00     5091     4410
## z_age_months        -0.38      1.77    -4.58     2.83 1.00     7242     2973
## method1             -0.41      3.54    -9.07     5.50 1.01     4918     1421
## IDS_pref             0.66      6.25    -7.54    14.21 1.00     2752      968
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   155.31     10.18   136.85   176.99 1.00     7865     5540
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_1.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + 
                      (1 | labid),
                      data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_null_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =312)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '3d251476a88e0247333a7c190c25e48a' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 6.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.64 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.494055 seconds (Warm-up)
## Chain 1:                0.475653 seconds (Sampling)
## Chain 1:                0.969708 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '3d251476a88e0247333a7c190c25e48a' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 5.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.581779 seconds (Warm-up)
## Chain 2:                0.51256 seconds (Sampling)
## Chain 2:                1.09434 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '3d251476a88e0247333a7c190c25e48a' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.71 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.520601 seconds (Warm-up)
## Chain 3:                0.405477 seconds (Sampling)
## Chain 3:                0.926078 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '3d251476a88e0247333a7c190c25e48a' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 8.7e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.87 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.553748 seconds (Warm-up)
## Chain 4:                0.495358 seconds (Sampling)
## Chain 4:                1.04911 seconds (Total)
## Chain 4:
summary(lmer.bf_null_1.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + (1 | labid) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.84      0.64     0.03     2.37 1.00     5601     3361
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           49.67      3.21    43.64    56.03 1.00     4019     3933
## CDI.z_age_months    22.91      5.87    14.36    34.56 1.00     4085     6110
## gender1              0.56     12.09   -15.13    29.92 1.00     4974     4212
## z_age_months        -0.38      1.82    -4.81     2.88 1.00     6882     2726
## method1             -0.29      3.16    -7.45     4.94 1.00     5668     2022
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   155.74     10.33   137.24   177.67 1.00     7565     4852
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_1.uk <- bridge_sampler(lmer.bf_full_1.uk, silent = T)
marloglik_null_1.uk <- bridge_sampler(lmer.bf_null_1.uk, silent = T)

Bayes factor

BF_full_1.uk <- bayes_factor(marloglik_full_1.uk, marloglik_null_1.uk)
BF_full_1.uk
## Estimated Bayes factor in favor of x1 over x2: 0.94008

4.8 UK model (interaction between IDS test age and IDS preference)

lmer.bf_full_2.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                      (1 | labid),
                     data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=546)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '5b897001df2f84a1996c21892ffd87d5' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.654481 seconds (Warm-up)
## Chain 1:                0.801309 seconds (Sampling)
## Chain 1:                1.45579 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '5b897001df2f84a1996c21892ffd87d5' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7.2e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.72 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.04675 seconds (Warm-up)
## Chain 2:                0.558659 seconds (Sampling)
## Chain 2:                1.60541 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '5b897001df2f84a1996c21892ffd87d5' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000416 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.16 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.10452 seconds (Warm-up)
## Chain 3:                0.972033 seconds (Sampling)
## Chain 3:                2.07656 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '5b897001df2f84a1996c21892ffd87d5' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.6 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.784626 seconds (Warm-up)
## Chain 4:                0.637108 seconds (Sampling)
## Chain 4:                1.42173 seconds (Total)
## Chain 4:
summary(lmer.bf_full_2.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 | labid) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.84      0.64     0.03     2.39 1.00     3969     2427
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    49.50      3.30    43.15    55.87 1.00     3525
## CDI.z_age_months             22.75      5.82    14.30    34.33 1.00     3826
## gender1                       1.05     12.46   -15.70    32.43 1.00     4139
## z_age_months                 -0.40      1.78    -4.79     2.79 1.00     6214
## method1                      -0.25      3.32    -8.11     5.76 1.01     6002
## IDS_pref                      0.25      4.67    -8.55    10.43 1.00     5497
## method1:IDS_pref             -0.16      4.74    -9.75     8.10 1.00     5314
## CDI.z_age_months:IDS_pref     0.32      3.38    -5.96     8.16 1.00     4842
## z_age_months:IDS_pref         0.10      3.29    -6.72     7.44 1.00     4534
##                           Tail_ESS
## Intercept                     4104
## CDI.z_age_months              6123
## gender1                       3799
## z_age_months                  3105
## method1                       1824
## IDS_pref                      1840
## method1:IDS_pref              1657
## CDI.z_age_months:IDS_pref     1953
## z_age_months:IDS_pref         1949
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   155.74     10.44   136.75   177.77 1.00     7335     5236
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_2.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + 
                      (1 | labid),
                      data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_null_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =689)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'fed304567041233af949a5c50d56c33a' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.45 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.776728 seconds (Warm-up)
## Chain 1:                0.942291 seconds (Sampling)
## Chain 1:                1.71902 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'fed304567041233af949a5c50d56c33a' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.811237 seconds (Warm-up)
## Chain 2:                0.883222 seconds (Sampling)
## Chain 2:                1.69446 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'fed304567041233af949a5c50d56c33a' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 5.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.53 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.744778 seconds (Warm-up)
## Chain 3:                0.957237 seconds (Sampling)
## Chain 3:                1.70202 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'fed304567041233af949a5c50d56c33a' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.63 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.18464 seconds (Warm-up)
## Chain 4:                1.00483 seconds (Sampling)
## Chain 4:                2.18947 seconds (Total)
## Chain 4:
summary(lmer.bf_null_2.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + (1 | labid) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.84      0.64     0.03     2.40 1.00     4572     2694
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    49.52      3.36    43.12    56.13 1.00     3236
## CDI.z_age_months             22.83      5.80    14.36    34.30 1.00     3651
## gender1                       0.84     12.49   -15.51    31.91 1.00     5130
## z_age_months                 -0.40      1.80    -4.84     2.78 1.00     6954
## method1                      -0.28      3.56    -8.65     6.31 1.00     5584
## IDS_pref                      0.52      6.32    -7.96    12.25 1.01     4157
## method1:IDS_pref             -0.10      5.19    -9.43     8.13 1.01     6863
## CDI.z_age_months:IDS_pref     0.24      2.86    -5.29     6.84 1.00     6084
##                           Tail_ESS
## Intercept                     3398
## CDI.z_age_months              5763
## gender1                       4107
## z_age_months                  3011
## method1                       1793
## IDS_pref                      1242
## method1:IDS_pref              2120
## CDI.z_age_months:IDS_pref     2366
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   155.84     10.39   136.89   177.51 1.00     8376     5386
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_2.uk <- bridge_sampler(lmer.bf_full_2.uk, silent = T)
marloglik_null_2.uk <- bridge_sampler(lmer.bf_null_2.uk, silent = T)

Bayes factor

BF_full_2.uk <- bayes_factor(marloglik_full_2.uk, marloglik_null_2.uk)
BF_full_2.uk
## Estimated Bayes factor in favor of x1 over x2: 0.91752

4.9 Combined model (interaction between lanaguage zone and IDS preference)

lmer.bf_full_3.combined <- brm(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + IDS_pref:language_zone +
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                           (1 + CDI.z_age_months | labid),
                         data = data.uk_nae,
                        family = gaussian,
                      prior = prior_combined_full_info,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .9999), 
                      seed=100)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '2ae4a1ae26a91de10302d79941e0a307' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000582 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.82 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 12.3233 seconds (Warm-up)
## Chain 1:                7.96945 seconds (Sampling)
## Chain 1:                20.2927 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '2ae4a1ae26a91de10302d79941e0a307' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000139 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.39 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 11.6286 seconds (Warm-up)
## Chain 2:                8.87576 seconds (Sampling)
## Chain 2:                20.5044 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '2ae4a1ae26a91de10302d79941e0a307' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000454 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.54 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 11.724 seconds (Warm-up)
## Chain 3:                8.00645 seconds (Sampling)
## Chain 3:                19.7304 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '2ae4a1ae26a91de10302d79941e0a307' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000173 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.73 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 12.0478 seconds (Warm-up)
## Chain 4:                13.6773 seconds (Sampling)
## Chain 4:                25.7251 seconds (Total)
## Chain 4:
summary(lmer.bf_full_3.combined)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months + method + IDS_pref + IDS_pref:language_zone + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 + CDI.z_age_months | labid) 
##    Data: data.uk_nae (Number of observations: 401) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 12) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.07      0.04     0.01     0.16 1.00
## sd(CDI.z_age_months)                0.03      0.01     0.02     0.06 1.00
## cor(Intercept,CDI.z_age_months)     0.14      0.38    -0.61     0.78 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       1628     1950
## sd(CDI.z_age_months)                3157     5279
## cor(Intercept,CDI.z_age_months)     1855     3392
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     0.33      0.03     0.28     0.39 1.00     4711
## CDI.z_age_months              0.04      0.01     0.02     0.06 1.00     3982
## language_zone1                0.09      0.04     0.02     0.16 1.00     3801
## gender1                       0.04      0.01     0.02     0.06 1.00     9997
## z_age_months                 -0.00      0.00    -0.01     0.01 1.00     8365
## method1                      -0.00      0.03    -0.07     0.07 1.00     4194
## method2                      -0.00      0.05    -0.12     0.09 1.00     3421
## IDS_pref                     -0.01      0.04    -0.08     0.07 1.00     8923
## language_zone1:IDS_pref       0.07      0.05    -0.04     0.17 1.00     7604
## method1:IDS_pref             -0.02      0.05    -0.12     0.08 1.00     7668
## method2:IDS_pref             -0.05      0.08    -0.21     0.11 1.00     6009
## CDI.z_age_months:IDS_pref     0.01      0.01    -0.01     0.03 1.00     9258
## z_age_months:IDS_pref        -0.00      0.01    -0.02     0.02 1.00     8977
##                           Tail_ESS
## Intercept                     3933
## CDI.z_age_months              4494
## language_zone1                3332
## gender1                       5754
## z_age_months                  6109
## method1                       3863
## method2                       4142
## IDS_pref                      6434
## language_zone1:IDS_pref       5949
## method1:IDS_pref              6349
## method2:IDS_pref              6147
## CDI.z_age_months:IDS_pref     6391
## z_age_months:IDS_pref         6077
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.21      0.01     0.20     0.23 1.00     8268     5679
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_3.combined <- brm(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + 
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                           (1 + CDI.z_age_months | labid),
                         data = data.uk_nae,
                      family = gaussian,
                      prior = prior_combined_null_info,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .9999), 
                      seed =200)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL '1e9fb42f69dff2e35aab04bc0c20c833' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000615 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.15 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 13.2027 seconds (Warm-up)
## Chain 1:                8.69344 seconds (Sampling)
## Chain 1:                21.8962 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '1e9fb42f69dff2e35aab04bc0c20c833' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000373 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.73 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 12.3354 seconds (Warm-up)
## Chain 2:                13.267 seconds (Sampling)
## Chain 2:                25.6023 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '1e9fb42f69dff2e35aab04bc0c20c833' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000189 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.89 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 12.4057 seconds (Warm-up)
## Chain 3:                7.8485 seconds (Sampling)
## Chain 3:                20.2542 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '1e9fb42f69dff2e35aab04bc0c20c833' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000322 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 3.22 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 11.9008 seconds (Warm-up)
## Chain 4:                8.74046 seconds (Sampling)
## Chain 4:                20.6412 seconds (Total)
## Chain 4:
summary(lmer.bf_null_3.combined)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 + CDI.z_age_months | labid) 
##    Data: data.uk_nae (Number of observations: 401) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 12) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.06      0.04     0.01     0.15 1.00
## sd(CDI.z_age_months)                0.03      0.01     0.02     0.06 1.00
## cor(Intercept,CDI.z_age_months)     0.11      0.39    -0.67     0.77 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       1755     3160
## sd(CDI.z_age_months)                3490     4973
## cor(Intercept,CDI.z_age_months)     1652     3129
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     0.33      0.03     0.28     0.39 1.00     4857
## CDI.z_age_months              0.04      0.01     0.02     0.06 1.00     3471
## language_zone1                0.09      0.03     0.02     0.16 1.00     3922
## gender1                       0.04      0.01     0.02     0.06 1.00    10010
## z_age_months                 -0.00      0.00    -0.01     0.01 1.00     8804
## method1                      -0.00      0.03    -0.07     0.06 1.00     4326
## method2                      -0.00      0.05    -0.11     0.09 1.00     3617
## IDS_pref                     -0.01      0.04    -0.08     0.07 1.00     7461
## method1:IDS_pref             -0.02      0.05    -0.12     0.09 1.00     6112
## method2:IDS_pref              0.01      0.06    -0.11     0.14 1.00     6422
## CDI.z_age_months:IDS_pref     0.01      0.01    -0.01     0.03 1.00     9622
## z_age_months:IDS_pref        -0.01      0.01    -0.03     0.02 1.00     7741
##                           Tail_ESS
## Intercept                     3521
## CDI.z_age_months              4281
## language_zone1                3649
## gender1                       6022
## z_age_months                  6358
## method1                       4570
## method2                       4296
## IDS_pref                      6270
## method1:IDS_pref              6007
## method2:IDS_pref              5734
## CDI.z_age_months:IDS_pref     5430
## z_age_months:IDS_pref         5676
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.21      0.01     0.20     0.23 1.00     8717     5608
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_3.combined <- bridge_sampler(lmer.bf_full_3.combined, silent = T)
marloglik_null_3.combined <- bridge_sampler(lmer.bf_null_3.combined, silent = T)

Bayes factor

BF_full_3.combined <- bayes_factor(marloglik_full_3.combined, marloglik_null_3.combined)
BF_full_3.combined
## Estimated Bayes factor in favor of x1 over x2: 0.09910

```